ls#!/bin/bash

ecut_energys=( 90 100 110 120 130 ) # the unit is Rydbergand 1 Rydberg = 13.6 eV

cat > ecut_energy_time << EOF
ecut    energy/atom    Time
EOF


echo "start to test ecut with energy and time"
for ecut_energy in ${ecut_energys[@]}
    do
        mkdir "ecut_${ecut_energy}"
        cd ecut_${ecut_energy}

        echo "start to write the INPUT of ecut ${ecut_energy} in path $PWD"
        cat > INPUT << EOF
INPUT_PARAMETERS
#Parameters (1.General)
calculation         scf
#Parameters (2.Iteration)
ecutwfc             100
scf_nmax            100
symmetry            0
kspacing            0.08
#Parameters (3.Basis)
basis_type          pw
#Parameters (4.Smearing)
smearing_method     gauss
smearing_sigma      0.002
#Parameters (5.Mixing)
mixing_type         pulay
mixing_beta         0.7
cal_stress          1
cal_force           1
EOF

echo "start to write the KPT of ecut ${ecut_energy} in path $PWD"
        cat > KPT << EOF
K_POINTS
0
Gamma
6 6 6 0 0 0
EOF

            cp /home/deep/DP-B-Si/DP-B-Si/confs/Si/STRU/Si8.STRU ./STRU
            cp /home/deep/DP-B-Si/DP-B-Si/potential势函数/Si_ONCV_PBE-1.0.upf ./Si_ONCV_PBE-1.0.upf

            echo "start to run abacus in path $PWD"
            OMP_NUM_THREADS=1 mpirun -n 5 abacus | tee out.log

            Time=`grep "TOTAL  TimeE" out.log | awk '{printf int($4)}'`
            energy_atom=`grep "FINAL_ETOT_IS" OUT.ABACUS/running_scf.log | awk -v stru="8" '{printf "%.8f\n",$2/stru}'`
            echo "$ecut_energy   $energy_atom   $Time" >>../ecut_energy_time
        cd ..
done

echo "start to write the python script in path $PWD"

cat > ecut_energy_time.py << EOF
import matplotlib.pyplot as plt

with open('ecut_energy_time', 'r') as f:
    data = f.readlines()
print("strat read the data")
ecut = []
energy_atom = []
Time = []

for line in data[1:]:
    line = line.strip().split()
    ecut.append(str(line[0]))
    energy_atom.append(float(line[1]))
    Time.append(int(line[2]))

fig, (ax1, ax2) = plt.subplots(2, 1, sharex=True)
ax1.plot(ecut, energy_atom, marker='o')
ax1.set_ylabel('energy/atom')
ax1.set_title('Plot of ecut vs energy/atom')
ax1.set_ylim(min(energy_atom), max(energy_atom))

ax2.plot(ecut, Time, marker='o')
ax2.set_xlabel('ecut')
ax2.set_ylabel('Time')
ax2.set_title('Plot of ecut vs Time')
ax2.set_ylim(min(Time), max(Time))
plt.savefig('ecut_energy_time.png')

EOF

python ecut_energy_time.py

echo "all finish"






